Análisis Espacial

Minicurso SOCHE

Autora

Stephanie Orellana Bello

1 Análisis de patrones de puntos

1.1 Cargar datos

Trabajaremos con un set de datos que descargaremos con el paquete {rgbif}.

GBIF —Infraestructura Mundial de Información en Biodiversidad— es una organización internacional y una red de datos financiada por gobiernos de todo el mundo, destinada a proporcionar a cualquier persona, en cualquier lugar, acceso abierto y gratuito a datos sobre cualquier tipo de forma de vida que hay en la Tierra.

ESto está en formato Darwin Core:

El set de datos también se encuentra disponible en el github.

# install.packages("rgbif")

data <- rgbif::occ_data(scientificName = "Jubaea chilensis", 
                        limit=10000, 
                        country = "CL", 
                        basisOfRecord = "HUMAN_OBSERVATION")$data

1.2 Convertir a objeto espacial

Al descargar este set de datos, estará disponible como un data.frame por lo que es necesario transformarlo a un objeto espacial para hacer algunas operaciones. Esto lo haremos con el paquete {sf}:

# install.packages("sf")
# install.packages("tidyverse")

library(sf)
library(tidyverse)

data_sp <- data %>% 
  filter(!is.na(decimalLongitude), !is.na(decimalLatitude)) %>% 
  st_as_sf(coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)

1.3 Visualizar

Podemos visualizarlo de manera rápida con el paquete {mapview}:

# install.packages("mapview")
library(mapview)

mapview(data_sp)

Ahora podemos comenzar con nuestro análisis de patrones de puntos:

1.4 Corroborar sistema de proyección

Este tipo de análisis solo puede ser hecho con sistemas de coordinadas proyectados, entonces es necesario cambiar nuestro sistema de coordenadas:

data_sp_utm <- data_sp %>% 
  st_transform(crs = 32719)
data_sp_utm$geometry
Geometry set for 360 features 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 221553.9 ymin: 6075129 xmax: 356415 ymax: 6692283
Projected CRS: WGS 84 / UTM zone 19S
First 5 geometries:
POINT (259482.3 6337580)
POINT (259505 6337556)
POINT (259517.3 6337518)
POINT (341457.1 6236730)
POINT (258941.6 6337451)
st_crs(data_sp_utm)
Coordinate Reference System:
  User input: EPSG:32719 
  wkt:
PROJCRS["WGS 84 / UTM zone 19S",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 19S",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-69,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",10000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Between 72°W and 66°W, southern hemisphere between 80°S and equator, onshore and offshore. Argentina. Bolivia. Brazil. Chile. Colombia. Peru."],
        BBOX[-80,-72,0,-66]],
    ID["EPSG",32719]]
st_crs(data_sp)
Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

1.5 Establecer ventana de análisis

Como primera aproximación utilizaré la extensión de los puntos:

win <- data_sp_utm %>% 
  st_bbox() %>% 
  st_as_sfc()

mapview(list(data_sp_utm, win))

1.6 Crear patrón de puntos

Ahora utilizaremos el paquete {spatstat} para realizar nuestro análisis. Para crear un patrón de puntos es necesario tener las coordenadas de los puntos y la ventana de análisis:

# install.packages("spatstat")

library(spatstat)

x <- st_coordinates(data_sp_utm)[,1]
y <- st_coordinates(data_sp_utm)[,2]
box <- as.owin(win)   
  
pp1 <- ppp(x = x, y = y, window = box)
plot(pp1)

1.7 Realizar kernel de densidad

Para calcular la intensidad de uso utilizaremos la función density:

kernel <- density.ppp(pp1)
plot(kernel)

library(stars)
Warning: package 'stars' was built under R version 4.2.1
Loading required package: abind
kernel_stars <- kernel %>%  st_as_stars() %>% st_set_crs(32719)

mapview(kernel_stars)

Podemos estandarizar el kernel haciendo un poco de álgebra:

kernel_stars$valor_01 <- (kernel_stars$v - min(kernel_stars$v))/(max(kernel_stars$v)- min(kernel_stars$v)) 

mapview(kernel_stars[2])

1.8 K de ripley

La función K de Ripley, K(r) es un método basado en la distancia que mide la aglomeración de un patrón de puntos espacial contando el número medio de vecinos que presenta cada punto dentro de un círculo de radio (r) determinado, en un determinado espacio, la función K compara el valor observado a una cierta distancia con el valor esperado a esa misma distancia; dado un proceso de Poisson homogéneo, también conocido como complete spatial randomness (CSR), es decir, que todos los puntos tienen la misma probabilidad de ocurrir en cualquier parte del área de estudio Fuente

Hacemos el análisis de K de Ripley con la función Kest

k_ripley <- Kest(pp1)
plot(k_ripley)

Para más información sobre la interpretación, VER